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We present an algorithm for solving the general relativistic initial value equations for a corotating 
polytropic star in quasicircular orbit with a nonspinning black hole. The algorithm is used to 
obtain initial data for cases where the black hole mass is 1, 3, and 10 times larger than the mass 
of the star. By analyzing sequences of constant baryon mass, constant black hole mass initial data 
sets and carefully monitoring the numerical error, we find innermost stable circular orbit (ISCO) 
configuration for these cases. While these quasiequilibrium, conformally flat sequences of initial 
data sets are not true solutions of the Einstein equations (each set, however, solves the full initial 
value problem), and thus, we do not expect the ISCO configurations found here to be completely 
consistent with the Einstein equations, they will be used as convenient starting points for future 
numerical evolutions of the full 3+1 Einstein equations. 
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I. INTRODUCTION 

The inspiraling black hole (BH) / neutron star (NS) 
binary system is thought to be both a promising candi- 
date for the central engines of gamma-ray bursts (GRBs) 
and a likely source of gravitational waves detectable by 
ground-based laser interferometric gravitational wave de- 
tectors (LIGO/VIRGO/TAMA/GEO) or space-based in- 
terferometers (LISA). Theoretical estimates of the event 



rate of BH/NS mergers [pj give 



10- 



per year per 



galaxy. For the advanced detectors of LIGO (sensitive 
out to 1000 Mpc), this would amount to approximately 
1 event per day. It is quite likely that LIGO-II will be 
able to extract information about the equation of state of 
the nuclear matter inside a neutron star from the gravi- 
tational waves of BH/NS mergers In addition, hard 
GRBs with relatively short duration [|| could be pro- 
duced from BH/NS mergers @-§. 

Apart from the fact that the BH/NS binary system is 
an interesting problem in and of itself, a better under- 
standing of the details of BH/NS mergers is important 
for observational general relativistic astrophysics. While 
some aspects of the BH/NS merger have been investi- 
gated with Newtonian simulations 0,^, to date, the com- 
plexity of both the Einstein equations and the relativistic 
hydrodynamical equations have prevented studies involv- 
ing full general relativistic simulations of BH/NS merg- 
ers. Here, we present the first steps in the direction of 
a fully general relativistic treatment of the coalescence 
of an inspiraling binary BH/NS, namely, the calculation 
of general relativistic initial data corresponding to a BH 
and a quasiequilibrium NS in a quasicircular orbit, near 
the innermost stable circular orbit (ISCO). In contrast 
to Newtonian theory, initial data in general relativity is 



not arbitrary in that one does not have complete freedom 
in the initial choices of the dynamical fields. The gravi- 
tational field (specified by the 3-metric) and the matter 
fields must obey the Hamiltonian and momentum con- 
straints that take the form of four coupled, elliptic, par- 
tial differential equations. We present an algorithm for 
numerically solving these constraint equations, together 
with conditions that specify initial data corresponding 
to a quasiequilibrium polytropic neutron star in a coro- 
tational orbit with a black hole. While such a treatment 
has been carried out for NS/NS (gj^ and BH/BH sys- 
tems |l^Jl^ , this is the first such treatment for a general 
relativistic BH/NS system. In addition, we find config- 
urations that correspond to an approximate innermost 
stable circular orbit for the BH/NS system by defining 
an effective binding energy, and locating the minimum 
of this binding energy for constant rest mass, constant 
BH mass initial data sequences. Although these are only 
approximations to the true ISCO configurations (since 
a time series of these constant rest mass, constant BH 
mass sequences are not actually solutions to the Einstein 
equations), they will provide a place to begin studies in 
a full 3D numerical general relativistic treatment. 

The remainder of this paper is outlined as follows. In 
Section the equations and all assumptions are pre- 
sented. In Section HI, the numerical algorithm for solv- 
ing the equations for the BH/NS system is presented. 
In Section |^ we present a method of testing the code 
to guarantee that we are solving all equations correctly. 
In Section ^ we define an effective binding energy and 
calculate this binding energy for sequences of constant 
rest mass, constant BH mass initial data sets. The data 
set that corresponds to a minimum (when it exists) of 
the binding energy is taken to be an approximate ISCO 



1 



configuration. We briefly summarize our results in Sec- 
tion 0. 



II. EQUATIONS 

The equations one must solve to specify initial data 
for general relativistic systems are the full 4D Einstein 
equations 



Gab — SirTab 



(1) 



(here, we are using units where the gravitational con- 
stant G and the speed of light c are set equal to 1), pro- 
jected into a direction normal to a spatial Cauchy sur- 
face (we assume, as always in numerical relativity, the 
entire spacetime to be globally hyperbolic). These are 
the Hamiltonian and momentum constraints, given re- 
spectively by 



and 



(2) 



(3) 



which are elliptic conditions on the initial 3-metric jab 
and extrinsic curvature Kab- Here, ^'^^i? is the scalar cur- 
vature of the 3-metric, K is the trace of the extrinsic cur- 
vature Kab, n"' is the future directed unit vector normal 
to the Cauchy surface, and V is the covariant derivative 
operator compatible with the 3-metric jab- While it is 
not necessary to specify the lapse and shift functions as 
part of the initial data, we find it convenient to do so 
here. We take the conditions on the lapse a and shift 
vector [}°- to be, respectively. 



and 



CtK = Q 



V'llab = 



where Y,ab is the distortion tensor, defined by 



1 



(4) 



(5) 



(6) 



Eq. i is the maximal slicing condition for the lapse, while 
Eq. 5 is the minimal distortion equation for the shift. 

We now make several simplifying assumptions on the 
form of the metric jab, extrinsic curvature Kab, and mat- 
ter variables that make up the matter stress energy tensor 
Tab- First, introduce Cartesian spatial coordinates {x*}, 
where the Latin indices vary from 1 to 3. We assume the 
spatial 3-metric to be conformally flat 



(7) 



and the trace of the extrinsic curvature to vanish 



K = Ki-iY^ = 



(8) 



For convenience, we introduce the conformal tracefree 
extrinsic curvature A'^ to be 



(9) 



Note that we raise and lower conformal quantities with 
the conformal (flat) metric 6ij, so that, e.g., Aij = (fP'Kij. 
Furthermore, we assume the stress energy tensor to be 
that of a perfect fluid, namely. 



T 



ah 



pohu'^u'' + Pg 



ah 



(10) 



where po is the rest mass density of the fluid, u° is the 
4- velocity of the fluid, P is the pressure of the fluid, g"'' 
is the 4-metric, and h is the relativistic specific enthalpy 
given by 



h = 1 



(11) 



where e is the specific internal energy density of the fiuid. 
For the rest of this paper, we assume the perfect fluid 
equation of state, along with the polytropic equation of 
state: 



P = (F - l)poe = kpo 



(12) 



where F is the adiabatic index of the fluid (all numerical 
results for this paper use F = 2) and k is the polytropic 
constant of the fluid. For convenience, we completely 
flx the units by setting the polytropic constant fc, along 
with G and c, equal to 1. We deflne the 3- velocity with 
respect to velocity of the coordinates {a;*}, and note that 
these components are related to the spatial components 
of the 4- velocity, u', by the following: 

= W{v' - 13' /a) (13) 

where W is the Lorentz factor, W = au^ = 
(l-7y«V")"'/'. 



A. General Relativistic Equations 

There are several ways to accommodate the presence of 
a black hole in the computational domain for initial data 
preparation, including conformal imaging with horizon 
excision We adopt a method given by Brandt and 
Briigmann [p^ , where the singular nature of the black 
hole is taken explicitly into account by rewriting the con- 
formal factor (/) in terms of the function x as 



BH 



2 If - xbh\ 



X, 



(14) 



where Mbh is the bare mass of the black hole puncture 
(see |l^) and xbh is the coordinate location of the black 



2 



hole puncture. The Hamihonian constraint, Eqj2[ is now 
written as 

drd'x + Ir'AjA^' + 2TTcl>^pohW^ -P) = Q. (15) 

o 

As with the conformal factor, we rewrite the conformal 
tracefree extrinsic curvature A*-' in terms of the functions 
Z^i as 



2a \ 3 



(22) 



which is the condition that the time derivative of the 
conformal metric vanish, namely, Cti'P^^lij) = (in this 
case, the vector potential W' assumes the role of the 
shift). However, Eqs. |l6| and ^ reduce to Eq. |2| as 
Mbh ^ 0. 



2\x ~ xbh\ 

-{6'^ - n'BH'<^BH)PBHnkBH 



(16) 



where Pbh linear momentum of the black hole as 

measured from infinity, and = [x^ — x^^^) /\x — xbh\- 
The momentum constraint, Eq. ^ is now written as 



djZ'^ - 87r</)^Vo/iW^ V* = 0. 



Making the ansatz 

,7 



r 

2-q 



(17) 



(18) 



giving the six functions Z^^ in terms of the three functions 
W'' , the three momentum constraints can be written as 
equations for the three functions W^: 



3 ■' 'x^ 



'16T:r](t)^°X''^PohW^v' = 0. 



(19) 



Using the Einstein equations, the maximal slicing con- 
dition, Eq. ^ can be written in terms of the function 
T] = a4> as 

^^^'r] + 77 (-'^(j)-^AijA'^ + 27r(t)^{2poh - SpohW^ - 5P) 



0. 



(20) 



The minimal distortion shift equation, Eq. ^ can be writ- 
ten 



B. Matter Equations 

In writing down the equations that govern the matter, 
we follow closely the formalism of where corotating 
binary neutron stars were solved in the quasiequilibrium 
approximation. First, we assume that the 4-velocity vec- 
tor field of the fluid m° is proportional to an approximate 
Killing vector field, and write the 4-velocity as 



d 



(23) 



(recall that the vector field representing the flow of time 
is t"" = an"' + 13°'), where VI is the (constant) angular 
velocity. The normalization condition on the 4-velocity, 
u°'Ua — —1, now becomes 



2 ? 7 

V = ^ijV V-' = 



(3 + n 



(90) 



(24) 



The relativistic Bernoulli equation, under the assumption 



that the EOS is of the form of Eq. 12, can be directly 
integrated to yield 



constEint 



(25) 



which can now be written as 

a2 _ ^4 |^(^. _ y^f + + + (^.)2 j ^ Cb, 

(26) 

where Cb is a constant. 



+ 2(j)-^A'^dja + 16TTa(j)'^pohW^v' 



(21) 



Notice that the ansatz for the form of the extrinsic cur- 
Iq , is not equivalent to that used 
Due to the presence of a BH 
with non-zero momentum p6[ , we 
cannot put the extrinsic curvature in the form (analogous 
to Eq. |l|) 



vature, Eqs. and 
in NS/NS studies g 
in the form of Eq. 04 



III. NUMERICAL ALGORITHM 

We will now describe the numerical algorithm used to 
solve the equations in the previous section. What we re- 
quire are simultaneous solutions of the 9 equations (the 
Hamiltonian constraint (1), the maximal slicing condi- 
tion (1), the momentum constraints (3), the minimal dis- 
tortion shift equations (3), and the relativistic Bernoulli 
equation (1)) for the 9 fields x, W^, rj, (3^, and po, which 
completely specify the configuration. There are 3 free pa- 
rameters that completely specify the configuration of a 
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neutron star and black hole in quasicircular orbit about 
each other. These will be the the ratio of the mass of 
the black hole to the mass of the neutron star {fibn), the 
separation between the black hole and the neutron star 
(parametrized by the parameter xa, which is defined be- 
low), and the maximum rest mass density {po)rnax of the 
neutron star. The idea is to iteratively solve the required 
equations with a stable iterative procedure in which the 
residual of each equation decreases to some desired tol- 
erance. To this end, we start by following ^ in defining 
new spatial coordinates {x*} which are just a rescaling 
of the original coordinates {x*}: 



cr 



(27) 



We transform all quantities (except for the matter fields) 
to these new coordinates {i'}. All numerical computa- 
tions are done in these hatted coordinates. We can easily 
transform any calculated quantity back to the unhatted 

coordinates, e.g., di ~ di/a, A'^^ = A'^/cr, Q — f^/cr, or 

Mbh = <jMbh- 

We choose for the black hole and the neutron star to 
be orbiting in the 2 = plane, with their respective cen- 
ters lying on the x-axis. The neutron star surface will 
intersect the i-axis at two different points. We define 
the intersection furthest from the origin of the coordi- 
nate system to be exactly 1, and denote the intersection 
closest to the origin as xa- We require the bare mass of 
the black hole Mbh to be given in terms of the ADM 
mass of the neutron star 



Mbh = Pbn (Mq) 



ADM 



(28) 



where {Mo)adm denotes the ADM mass of a neutron star 
in isolation with total rest mass Mg. We further require 
the center of mass of the system to be at the origin of 
the coordinate system. We therefore set the coordinate 
of the black hole puncture Xgj^ — {xbh, 0, 0) to be 



along the x-axis. This is particularly easy in our case, 
as the conformal factor is written explicitly in terms of 
a I/tbh piece [tbh = the coordinate dis- 

tance from the black hole puncture) and X: which will 
be numerically solved with a robin boundary condition 
of(x-l) ~ 1/fjvs, where fjvs = - (xc, 0,0)* I . One 
could alternatively define the ADM mass of the neutron 
star {Mq)adm i'^ ^'^^ H terms of the l/fjv5 falloff of 
the function x, and one could argue that this would be 
a more accurate way of restricting the center of mass of 
the system to be at the origin of the coordinate system. 
We have checked that both methods of defining the ADM 
mass of the neutron star are the same to within 6% at 
worst, and much better than this in most cases. 



We determine the momentum of the black hole, P 



BH 



(0, P%fj, 0), from the condition that the total linear mo- 
mentum in the y-direction vanish. The linear momentum 
in the y-direction is given (see, e.g., EtI) by 



py 



Stt 



lim / (fSi {Ry' - ksy' 



(30) 



For the case of conformally flat initial data, where the 
conformal tracefree extrinsic curvature is written as in 
Eq. the condition that the total linear momentum in 
the y-direction vanish is simply 



pV 
^BH 



d?x a^(j)^°pohW^vy, 



(31) 



where the momentum constraint Eq. has been used to 
simplify the integral. 

The form of the Hamiltonian constraint (Eqs. p|,p^, 
momentum constraints (Eqs. p|Jl9|), and maximal slicing 
condition (Eqs. ^,^0|) that we solve numerically can now 
be written, respectively, as 

^^^'x + Icfy^'^A^jA^ + 2ncr^(f>^{pohW^ -P) = 0, (32) 



Xbh 



xc 

Pbn 



(29) 



where xc is the coordinate value along the x-axis where 
the rest mass density po of the neutron star acquires its 
maximum value. As xc will generally not be located 
on a discrete grid point, a quadratic polynomial is fit 
to the maximum discrete value and its nearest discrete 
neighbors in the x-direction. It is the coordinate loca- 
tion of the maximum of this polynomial that is denoted 
as Xc- One may worry that the prescription for the lo- 
cation of the black hole along the x-axis, Eq. |2^, may 
not actually place the center of mass of the system pre- 
cisely at the origin of the coordinate system. Analyti- 
cally, the center of mass of the system can only be de- 
termined by observers at infinity. In our case, one only 
needs to look at the 1 /f falloff of the conformal factor 



3 



and 



d,'^ - l6TTa^7j(l)^"x'''pohW^v' = 0, (33) 



d^d'r] + 77 ( -l<f)-^A^jA^ + 2Tra'^(f,'^{2poh - 3pahW^ - 5P) 



= 0. (34) 
In hatted coordinates, the conformal factor is written as 

Mbh 



2 \x-xbh\ 



X, 



(35) 



whereas the momentum constraint ansatz, Eq. |l8|, is writ- 
ten as 
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2r] \ 3 



(36) 



In numerically solving the matter equations (the rel- 
ativistic Bernoulli equation, Eq. p6| ), written in hatted 
coordinates here as 

h^a^ - 04 ((/3- - yClf + ipy + xClf + (/?-)') = Cb, 

(37) 

we again follow j|] by choosing three points along the x- 
axis, and evaluating Eq. ^ at these three points. Using 
Newton's method, we solve these three equations for the 
three constants cr, fi, and Cb- Using these three newly 
computed constants, we reset the matter variable po ev- 
erywhere through Eq. This uniquely specifies all of 
the matter variables h, v^, and W. The three points we 
evaluate Eq. |3^ are where the surface of the neutron star 
intersects the x-axis (x — 1,xa), and where po attains 
its maximum value (x = xc ). Solving Eq. || at three 
points for the constants cr, 51, and Cb is not straightfor- 
ward. While Cl and Cb appear explicitly in Eq. |3^, cr does 
not. However, the lapse a and the conformal factor (p de- 
pend on cr, as can be seen from Eqs. |3^ and In both 
of these equations, increasing (decreasing) a acts like an 
increase (decrease) in the matter source, thus resulting 
in a decrease (increase) in the lapse a, and an increase 
(decrease) in the conformal factor (/) in the vicinity of 
the NS. In the functions a — a{a) and (j) = 0(cr) 
were modeled using Newtonian scaling relations, which 
we found in our case to not produce a stable (converging) 
numerical algorithm for iteratively solving the differen- 
tial spacetime equations with the algebraic relativistic 
Bernoulli's equation. Instead, we use static spherical so- 
lutions of the Einstein's equations coupled to a perfect 
fluid (Tolman-Oppenheimer-Volkoff [|l|,0, or TOV, so- 
lutions) to model the functions a = a{a) and 4> — </*("■)■ 
For example, at i = xc, we take ac(cr), the function 
we will use to model the dependence of the lapse a on 
the scaling factor cr at the point of maximum rest mass 
density of the neutron star, to be 



ac{cr) = q;o + aToy(o') 



(38) 



where aTov{<^) is the lapse at the center of the TOV 
star with scaling factor cr, and ao is a constant specified 
by the consistency condition that the actual value of the 
lapse at X — Xc should be given by the original value of 
the scaling factor, aoid. 



ao = ac{(Joid) - aTov{<^oid) 



(39) 



Similar functions arc used to model the a dependence 
of the lapse function a at the surface of the stars {x — 
XA, 1), as well as for the conformal factor cf). We can now 
solve Eq. at the points x — xa,xcA using Newton's 



method, where the modeling functions (e.g., Eq. 38) are 
substituted for a and 0. 

The last equation we need to solve numerically is the 
minimal distortion shift equation, Eq. written in hat- 
ted coordinates as 
1 

r 



+ 2(f)'^A dja + 16TTa^a<t)^pohW^v\ (40) 

Numerically, we find the first term on the right hand side 
of Eq. ^ (terms involving {d0){d(j))/(t>) quite difficult to 
handle, due to the built in singular point of the conformal 
factor (j) at the black hole puncture, rBH = 0. However, 
if one writes down the minimal distortion shift equations 
in spherical coordinates about the point x^bh^ finds 
that the only nontrivial (non-constant) radial solutions 
to the equations are of the form /?* ~ rBH or /?* ^ 
as the radial coordinate about i^^, "Tbh = \x — xbh\i 
goes to zero. We therefore excise the region about the 
black hole for the minimal distortion shift equation. Note 
that, as the only place the shift appears (other than the 
minimal distortion shift equation) is in Bernoulli's equa- 
tion, Eq. H, which only applies where the matter content 
of the spacetime is nonzero, the minimal distortion shift 
equation is the only equation where black hole excision 
is needed. We excise a cubical region of the computa- 
tional domain for the minimal distortion shift equation, 
where the excision cube is centered about x^bh- '^^^ 
sides of the cube are taken to have a length of Mbh/'^ or 
3 Ai:, whichever is larger. We find that we need the cu- 
bical excised region to be at least three grid points wide 
in order to put a reasonable boundary condition there. 
The boundary condition that we use for the shift at the 
boundary of the excised region is 



/3* = -f 0{fBH) 



(41) 



where /Sg is a constant vector obtained by performing a 
Lorentz boost at infinity, with a boost factor correspond- 
ing to a velocity of Pbh / Mbh, on the Schwarzschild 
solution written in isotropic coordinates, as tbh ~^ 0- 
Specifically, we have 



/3S - 0, - 



Pbh/ Mbh 



1 - (Pbh/Mbh) 



= ,0 



(42) 



The algorithm is initialized by choosing values for the 
three parameters (po)maa;' ^^bm and xa, which will re- 
main fixed during the entire algorithm. We have coded a 
Full Approximation Scheme (FAS) multigrid solver |2^] 
to solve Eqs. |2|, and |4^ simultaneously, thus ob- 

taining solutions to the Hamiltonian constraint, momen- 
tum constraints, maximal slicing condition, and mini- 
mal distortion shift equations, respectively. We then 
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solve Bernoulli's equation, Eq. at the three points 
X = l,XA,xc for the three constants a, Cl, and Cb- 
These 3 constants are then used to set the matter field 
po using Bernoulli's equation Eq. We then set the 
constants Mbh, Pbh, and xbh via. the prescription 
Eqs. [ 



pl| , and |29[ respectively. We repeat the above 
process until the entire configuration has converged to a 
simultaneous solution. Typically, we stop the iteration 
process when the constants cr, fi, and Cb change by less 
than one part in 10,000 from one iteration to the next. 



IV. CODE VALIDATION 

An important and often overlooked ingredient in nu- 
merical relativity is code validation (e.g., see |21| and 
references therein). Due to the complexity of the Ein- 
stein equations, it is extremely important to verify that 
the coded finite difference equations are consistent with 
the differential equations one wants to solve. That is, as 
the discretization parameter A goes to zero, the finite 
difference equations that are solved in the code (denoted 
by the operator equation Ca = 0) should be verified to 
approach the differential equations (denoted by the op- 
erator equation £ = 0) as some integral power of A: 



(43) 



where the integer n is determined by the type of finite 
differencing used. Typically, one chooses second-order 
finite differencing (as we have done in this paper), in 
which case n = 2. Only when a numerical code used 
to solve Ca is shown to satisfy Eq. ^ can that code be 
considered validated. 

It turns out to be quite straightforward to verify any 
numerical solver Ca for Eq. using the notion of an 
Independent Residual Evaluator (IRE) ||2^. The idea is 
to construct, through an independent means (by either 
writing the differential equation in a different form or 
using different finite differencing schemes, or preferably 
both), a residual evaluator for the differential equations 
represented by C This is usually much easier to do than 
to code up a numerical solver for C. Denote this IRE of 
C to be C'a, and assume that the truncation error is of 
the same order n as that of the numerical solver of Ca 
implemented by the code. We now write Eq. W3 as 



Ca-C'a^O{A'^), 



(44) 



We can now compute numerically the left hand side of 
Eq. ^ directly. Simply use the code to solve the finite 
difference equations Ca = and compute the IRE C'a on 
the resulting numerical solution. By Eq. 44 , the value of 
this IRE should approach zero as ©(A"). By performing 
this numerical calculation at different resolutions (i.e., 
different values of the discretization parameter A), one 
checks whether or not this is the case. Only once the 



IRE is shown to approach zero as ©(A") can the code 
be considered to be validated. 

Of course, due to the specific nature of any particu- 
lar numerical solution, errors in a numerical code which 
solves the finite difference equations Ca may not become 
apparent with simply one particular numerical solution. 
This is why we advocate the verification of Eq. ^ not 
only with some fiducial numerical solution, but for all 
numerical results produced by a code\ 

In our case, we use a second order accurate Full Ap- 
proximation Scheme (FAS) multigrid method to solve the 
eight coupled elliptic partial differential equations which 
are the Hamiltonian constraint (Eq. |3^ ), the momen- 
tum constraints (Eqs. ^3| ), the maximal slicing condi- 
tion (Eq. |3^ ), and the minimal distortion shift equations 
(Eqs. 40). In addition, we code an IRE for each equation 
that we solve numerically. To guarantee that the IRE is 
truly independent of the numerical solver, we use the co- 
variant forms of the Hamiltonian constraint (Eq. mo- 
mentum constraints (Eqs. |^), maximal slicing condition 
(Eq. ^), and minimal distortion shift equations (Eqs. ||), 
computed in the hatted (computational) coordinates Sf 
with absolutely no assumptions on the form of the met- 
ric or extrinsic curvature. For each numerical solution 
computed using the process described in Section III, we 



calculate the value of each IRE, and verify that it goes to 
zero as 0{A?) (since we have used second-order accurate 
methods in both our solvers and our IREs). 

It would not be feasible to show the results of all of the 
convergence tests of the IREs which were, in fact, per- 
formed on each numerical solution obtained in compiling 
the results of this paper. Instead I will present an exam- 
ple of just one of the convergence tests of the IREs, and 
state that this same test was performed on all numeri- 
cal solutions obtained for this paper. For this example, 
we set the three free parameters as follows: the ratio of 
the mass of the black hole to the mass of the neutron 
star = 3, the separation parameter xa = 0.4668, 
and the maximum rest mass density of the neutron star 
(po),„Qa; — 0.09265. Using the algorithm described in 



Section III, we numerically solve for the 9 fields x, W^, 
rj, and po- This corresponds to a configuration in 
which the rest mass of the neutron star Mq is approx- 
imately 73% that of the maximum rest mass of a sta- 
ble TOV configuration (we have set the adiabatic index 
r = 2 for this example). Using this numerically gener- 
ated solution, we construct the 3-metric, extrinsic cur- 
vature, lapse, and shift. We then compute the IRE of 
the 8 partial differential equations at each point on the 
computational grid. This entire process is repeated for 
resolutions of Ax = 0.06, 0.03, and 0.015, using the num- 
ber of grid points rix along the i-axis to be = 64, 128, 
and 256, respectively. Shown in Figures gaud I are the 
values of the IRE along the i-axis at the three resolutions 
of the y-component of the momentum constraint, the y- 
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FIG. 1. An independent residual evaluator (IRE) evalu- 
ated on the numerical solution obtained for the configuration 
specified by fiBH = 3, = 0.4668, and {po)max ~ 0.09265. 
Shown is the value of the independent residual for the 
y-component of both the momentum constraint and min- 
imal distortion shift equation MDS^ along the i-axis for 3 
separate resolutions. 



FIG. 2. An independent residual evaluator (IRE) evalu- 
ated on the numerical solution obtained for the configuration 
specified by /ibh = 3, xa = 0.4668, and (po)^^^ = 0.09265. 
Shown is the value of the independent residual for the Hamil- 
tonian constraint H and maximal slicing equation K along 
the i-axis for 3 separate resolutions. 



component of the minimal distortion shift equation, the 
Hamiltonian constraint, and the maximal slicing condi- 
tion. Since we require the values of the IRE to approach 
zero as 0{Ax'^), we multiply the medium resolution IRE 
{Ax — 0.03) by a factor of 4, and we multiply the high 
resolution IRE {Ax — 0.015) by a factor of 16. In this 
way, it is simple to see if our numerical solver is consis- 
tent with the IRE (and thus consistent with the origi- 
nal differential equations): simply plot the values of the 
IREs scaled with the appropriate number (1,4, and 16 for 
the low, medium, and high resolutions, respectively) and 
see if the values are the same. If they are the same, we 
are assured that finite difference equations we are solving 
are consistent with the differential equations we want to 
solve, thus validating the code as described above. 

In observing the values of the IREs, we note that the 
residual appears not to be converging to zero as the sec- 
ond power of the discretization parameter Ax near the 
black hole, Xbh = —0.244. This is to be expected, as the 
conformal factor has a 1/f pole at that point, and higher 
order terms in the truncated Taylor series for the finite 
difference approximations are of the same order (or larger 
than) the second-order term near this pole. It could also 
be the case that higher order terms in the truncated Tay- 
lor series for the finite difference approximations of the 
y-component of both the momentum constraints and the 
minimal distortion shift equations (Figure 0) are caus- 
ing the lack of strict second-order convergence observed 
at the surface of the neutron star {x = 0.4668 and 1.0); 
both equations have source terms with a discontinuous 
derivative. It also could be the case that setting stronger 
tolerances for the multigrid elliptic solver would elimi- 



nate this behavior at these resolutions; we would nor- 
mally terminate the multigrid elliptic solve W-cycle pc| ] 
when the L2-norm of the residual fell below 5,000 times 
the L2-iiorm of the truncation error. 

Other than these noted deviations, we see strict 
second-order convergence of the IREs. This makes us 
confident that we have eliminated all errors in the code 
that could possibly cause the solutions to the difference 
equations to converge to anything other than solutions 
to the differential equations. 

For completeness, we also show an IRE for Bernoulli's 
equation, Eq in Figures || and ^. Although it is 
not a differential equation, and thus we do not expect 
to see convergence, we see that the algorithm described 
in Section [II is solving the relativistic Bernoulli's equa- 
tion to almost 1 part in 10,000 inside the neutron star, 
0.4668 <x < 1.0. 



V. SEQUENCES OF CONSTANT NS BARYONIC 
MASS, CONSTANT BH BARE MASS BH/NS 
BINARIES 

Now that we have successfully validated a code that 
can produce fully general relativistic initial data that cor- 
responds to a neutron star and a black hole in quasiequi- 
librium, quasicircular orbit, we are ready to use this as 
initial data for a numerical evolution code that can han- 
dle both black holes and general relativistic hydrodynam- 
ics. However, if we choose to numerically evolve an initial 
data configuration that is too far away from the inner- 
most stable circular orbit (ISCO), we will have to numer- 
ically follow many orbital periods before the final plunge. 
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FIG. 3. An independent residual evaluator (IRE) evalu- 
ated on the numerical solution obtained for the configuration 
specified by ^bh = 3, = 0.4668, and (po)^^^ = 0.09265. 
Shown is the independent residual of Bernoulli's equation, 
Eq. along the x-ajcis for 3 separate resolutions. We only 
require this equation to be satisfied where po 7^ 0, namely, for 
0.4668 < X < 1.0 (see Figure |). 
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FIG. 4. An independent residual evaluator (IRE) evalu- 
ated on the numerical solution obtained for the configuration 
specified by ^bh = 2>, xa = 0.4668, and (po)^^^. = 0.09265. 
Shown is the independent residual of Bernoulli's equation, 
Eq. along the i-axis for 3 separate resolutions. Only the 
X values where po 7^ are shown (0.4668 < i < 1.0). 



which would, at best, be a difficult task. Conversely, if 
we choose to numerically evolve an initial data configu- 
ration that is too far inside the ISCO, we will be starting 
with non-physical initial data. As stated in the abstract 
of this paper, the only way to truly find the ISCO of a 
particular system will be through fully general relativis- 
tic numerical studies. However, we can hope to find an 
approximate location of the ISCO by following |p|-|l3|| and 
studying sequences of initial data sets that have both con- 
stant baryonic mass and constant black hole mass. We 
will define an effective binding energy per unit rest mass, 
Eb, in terms of the ADM mass of the NS (Mjvs), the 
bare mass of the BH (Mbh), and the total ADM mass of 
the system Madm- We will take the minimum (when it 
exists) of this effective binding energy Ei, to be the con- 
figuration which approximates the ISCO configuration of 
the system. Fixing the rest mass of the system and the 
binding energy per unit mass in this way leaves one free 
parameter left from the initial 3 free parameters (i-ibm 
XA, and {po)max)- We can therefore find ISCO configu- 
rations for varying values of the ratio of the BH and NS 
mass, fj,b„. 

Obviously, the baryon number in the NS will be con- 
served during a quasiequilibrium orbit of a NS and BH. 
It may be argued, however, that one should fix a differ- 
ent mass (e.g., the apparent horizon mass) for the black 
hole when looking at sequences to try to find the ISCO. 
Technically, there is no exactly conserved BH mass dur- 
ing the quasiequilibrium orbit of a NS and BH; even the 
event horizon mass will not be strictly conserved during 
the quasiequilibrium orbits down to the ISCO. Here, we 
find that holding the bare mass of the BH fixed to be 
equivalent (to within numerical accuracy) to holding the 
apparent horizon mass fixed. For small values of fibn, this 
is due to the fact that the NS is very far away from the 
BH compared to the mass of the BH and thus does not 
much affect the area of the apparent horizon. For large 
values of fibn, while the NS is closer to the BH, the mass 
of the NS does not much affect the area of the apparent 
horizon. 

We define the binding energy per unit rest mass as 



Eb 



Madm - M 



BH 



NS 



(Mo, 



(45) 



NS 



where Mbh is the bare mass of the BH, (A/o)^5 is the 
rest mass of the NS, Mjvs is the ADM mass of a NS in 
isolation with a rest mass of (Mo)jy5, ^^'^ Madm is the 
ADM mass of the entire system, defined to be 



Madm = —r- lim 



1 

IGtt 



3 

1,3 = 1 



33 



V dx^ dx^ 



(46) 



where the integration surface is of constant radial coordi- 
nate r. Using Stoke's theorem, along with our particular 
form of the conformal factor (Eq. ^) of our conformally 
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flat 3-metric, we can write this integral for the ADM mass 
of the system as 



Madm = Mbh [ d^x ^'^^x 

Ztt J 



(47) 



The integrand d'^diX can be replaced using the Hamilto- 
nian constraint, Eq. ^ in terms of the conformal extrin- 
sic curvature, conformal factor, and matter source terms. 
It is in this form that we numerically integrate (in hatted 
coordinates a;*) to obtain the ADM mass of the system. 

Recall that, for any given mass ratio parameter /if,„, 
our numerical algorithm described in Section [II does 
not allow us to fix the rest mass of the NS. Instead, the 
remaining freedom of the system is parameterized in the 
numerical algorithm by the separation parameter xa and 
the maximum rest mass density (po)maa: of the NS. In or- 
der to construct sequences of constant rest mass initial 
data for a particular mass ratio /ib„, we must compute 
initial data sets for various values of parameters xa and 
iPo)max- interpolating these results, we can find con- 
stant rest mass sequences. Specifically, we chose various 
values (five, typically) for the separation parameter xa- 
Then, for each choice of xa we solve for the initial data 
set corresponding to a series of different parameter val- 
ues of {po)max ^^"i connect these results using a cubic 
spline (whose independent parameter is the rest mass of 
the system). We can then, for each value of xa, find the 
details of any initial data set determined by any partic- 
ular value of the rest mass. In this way, we can plot the 
binding energy per unit rest mass Ef, as a function of 
angular frequency parameter 

In Figure]^ we plot, for mass ratio parameter /x;,„ = 1, 
the binding energy per unit mass Ei, for a series of con- 
stant rest mass sequences. We consider the minimum 
of each sequence to be approximations to the ISCO 
configuration for the binary BH/NS system with the 
given rest mass. Note that no minimum exists for the 
Mo = 0.45(Mo)^^^ sequence in Figure ||. We do not 
claim that a ISCO configuration does not exist for a 
BH/NS system where Mq = 0.45(Mo)„^^. Instead, this 
is a result of the fact that our numerical algorithm de- 
scribed in Section HI does not converge on a solution 
when Xa < (i.e., when the center of mass of the system 
is inside the NS). Note from Table | that the smallest 
value of XA used in the construction of Figure || is, in 
fact, XA = 0. These points correspond to the highest 
value of O for each constant rest mass sequence in Fig- 
ure H These points are such that the center of mass of 
the system exists at the surface of the NS. 

Observe the error bars in Figure ||. It seems that every 
branch of science except numerical relativity has, in the 
past, at least tried to quantify the errors in any given 
measurement and/or calculation. While in the past it 
may have been true that any reasonable measure of er- 
ror in a 3D numerical relativity calculation would result 
in error bars of the same magnitude of (or larger than) 



-0.004 
-0.006 




FIG. 5. The binding energy per unit rest mass vs. the 
angular velocity (scaled by the rest mass of the system) for 
constant rest mass sequences where the ratio of the BH mass 
to the NS mass, ^ibn, is 1. Shown are sequences whose rest 
mass are 0.45, 0.6, 0.75, and 0.9 times the maximum stable 
rest mass of a NS for a F = 2 polytropic equation of state. 
In units where G — c = k = 1, the maximum rest mass 
for a r = 2 polytropic star is (A/q)^^^ = 0.179862. The 
minimum of each sequence (when it exists) is taken to be an 
approximate ISCO configuration (see Tables |^ and |l|). 





Tlx 


Ax 


XA 


1 


256 


0.015 


0.0, 0.08397, 


1 


128 


0.03 


0.1881, 0.3215, 


1 


160 


0.015 


0.5 


3 


256 


0.015 


0.3, 0.3771, 


3 


128 


0.03 


0.4668, 0.5727, 


3 


160 


0.015 


0.7 


10 


256 


0.01 


0.5, 0.5796, 


10 


128 


0.02 


0.6707, 0.7761, 


10 


160 


0.014375 


0.8358, 0.9 



TABLE I. Configurations used to estimate the truncation 
errors and boundary errors of the constant rest mass, constant 
BH mass sequences in Figures |1 |, and and of the ISCO 
configurations of Tables O and n| (see discussion) . 
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the actual quantity being calculated, this is no longer the 
case. Available computer size and power are now such 
that we can make reasonable attempts to quantify the er- 
rors in our numerical relativity calculations. In our case, 
there are two sources of numerical approximation errors. 
The first source of numerical error is usually referred to 
as 'truncation error', and is the error induced by neglect- 
ing higher order terms in the Taylor series expansion of 
our finite difference approximations of derivatives of func- 
tions. In our case, we have used second order differenc- 
ing, which means the highest order non-zero truncation 
error is of order Ai^, i.e., the truncation error goes to 
zero as Ai^ goes to zero. The second source of numeri- 
cal error is due to the fact that our computational outer 
boundary (which is cubical) is at a finite distance from 
the origin, whereas the boundary conditions to our ellip- 
tic equations are formulated in terms of the behavior of 
the fields at oo (e.g., (0 — 1) — > 0{l/f) as f — > oo). We 
use robin boundary conditions in our multigrid elliptic 
solver. E.g., we could place the condition (</> — 1) ^ (1/?^) 
on the field at a finite distance from the origin of our 
coordinate system. In doing so, we neglect higher order 
terms in a 1/f expansion of cj) at the boundary. Thus, the 
error we make by placing our computational boundary at 
a finite distance fc from the origin is of order ©(l/f^). 
Therefore, the total error made in the computation of 
any quantity Q in our numerical method can be modeled 
as the following relationship between the exact value Qe 
(the exact value one would obtain if one solved the differ- 
ential equations directly) and the quantity obtained by 
numerical calculation Qn as 



Qn — 



1 



(48) 



where et and Cb are constants and "• • •" represents the 
higher order contributions to the error. By performing 
a numerical calculation with 3 different sets of resolu- 
tion/gridpoints, we can solve for the unknown quantities 
Qe, et, and Cf,. We now have a measure of the numerical 
error for our best calculation Q„, namely, the calculation 
with the highest resolution and/or with computational 
boundary furthest from the origin. For the numerical 
calculations performed for Figure I (see Table | for the 
numerical configurations used in the calculation of the er- 
rors), as well as for all of the numerical results presented 
in this paper, we attach errors that are equal to 



(Ag), 



max{\Qr, 



1 



,|,|etAx^|,|efc-:^|} (49) 



where Qn, Ax, and fc is the result and configuration of 
the highest resolution run (the configuration with — 
256 in Table ^; Q^, et, and Cf, are computed (separately 
for each quantity Q) as described above. Note that, while 
this prescription (Eq. p9| ) for the size of the error bars 
can overestimate the error, it will never underestimate 
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FIG. 6. The binding energy per unit rest mass vs. the 
angular velocity (scaled by the rest mass of the system) for 
constant rest mass sequences where the ratio of the BH mass 
to the NS mass, /ii,„, is 3. Shown are sequences whose rest 
mass are 0.45, 0.6, 0.75, and 0.85 times the maximum stable 
rest mass of a NS for a F = 2 polytropic equation of state. 
In units where G — c — k = 1, the maximum rest mass 
for a F = 2 polytropic star is (Mo)„^^ = 0.179862. The 
minimum of each sequence (when it exists) is taken to be an 
approximate ISCO configuration (see Tables and |l]). 



the magnitude of the largest error terms in either the 
truncation error or the boundary error. 

In Figures |and@, we show the binding energy per unit 
mass Eb for a series of constant rest mass sequences for 
mass ratio ^bn = 3 and 10, respectively. Note from Ta- 
ble H that the minimum separation parameter xa used for 
the fibn = 3 and /ifc„ — 10 calculations are xa — 0.3 and 
XA — 0.5, respectively. We have found that, for lower val- 
ues of XA (i.e., smaller separation), the Roche limit |2^] 
is reached for large values of ipo)max- This result is prob- 
ably highly dependent on the EOS (we use an adiabatic 
index F = 2 here). However, as can be seen from Fig- 
ures ^ ^, and 0, a minimum in binding energy per unit 
rest mass Eb is found before the Roche limit is reached, 
at least for higher rest mass cases. Notice also that the 
maximum rest mass reported for Figures ^, ^, and |^ 
(corresponding to = 1,3, and 10, respectively) is 
0-9(A^o)™,,, 0.85(Mo)„,„ and 0.8(Mo)_„ respectively. 
We have found that, in each case, for slightly higher rest 
mass configurations, the numerical algorithm described 



in Section [II fails to converge to a simultaneous solution 
in a reasonable amount of iterations. However, since we 
are really only interested in initial data for future numer- 
ical simulations, a maximum rest mass of the NS in the 
range of 0-HMo)max " ^■^i^o)rnax certainly adequate. 

Note that the errors for the binding energy per unit 
mass Eb for the fibn = 10 case (Figure 0) are over an 
order of magnitude larger than the errors in the fXbn = 1 
case (Figure H). This is due to the fact that the relative 
separation is much larger in the higher fibn case (the xa 
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M„ = 0.45 (M„), 



FIG. 7. The binding energy per unit rest mass vs. the 
angular velocity (scaled by the rest mass of the system) for 
constant rest mass sequences where the ratio of the BH mass 
to the NS mass, Hbn, is 10. Shown are sequences whose rest 
mass are 0.45, 0.6, and 0.8 times the maximum stable rest 
mass of a NS for a F = 2 polytropic equation of state. In units 
where G = c= fc = l,the maximum rest mass for a F = 2 
polytropic star is (Mo)^^^ = 0.179862. The minimum of each 
sequence (when it exists) is taken to be an approximate ISCO 
configuration (see Tables H and III). 



used for /ift„ = 10 is closer to 1, see Table |); more reso- 
lution would be needed to get the same number of points 
across the NS. 

In Tables || and flU, we tabulate the results of the 



approximate ISCO configurations which we define to be 
the minimum of the binding energy per unit mass Eb 
for a constant rest mass sequence (Figures ^ - 0). To 
find the minimum, we fit a second order polynomial to 
the minimum data point, e.g. in Figure|5[ and its nearest 
neighbor on cither side. The minimum of this polynomial 
designates a good approximation of the minimum of the 
actual curve. The errors in Tables |l| and |!| 



are com- 



puted using the same method as that used to calculate 
the error bars in Figures || - (see Eq. . Note that the 
relative errors of the tabulated values increase for larger 
values of the BH/NS mass ratio parameter fitn- This is 
due to the fact that as the mass of the BH increase rel- 
ative to the mass of the NS, the separation between the 
two bodies for the ISCO configuration increases in units 
of the NS radius. As a result, there is less numerical reso- 
lution per NS radius for higher /if,„, resulting in a relative 
loss of accuracy. For configurations where fibn > 10, it 
is likely that mesh refinement techniques will be need to 
accurately resolve both the NS and BH. 





Mo/(Mo)_ 


Eb 


J /{Mo') 


nMo 




1 


0.6 


-0.00283 ± 0.00012 


2.9069 ± 0.0033 


0.00855 ± 0.00020 


0.06326 ± 0.00018 


1 


0.75 


-0.00353 ± 0.00025 


2.609 ±0.025 


0.01127 ±0.00099 


0.09589 ±0.00061 


1 


0.9 


-0.00403 ± 0.00024 


2.3887 ± 0.0039 


0.01417 ± 0.00038 


0.14947 ±0.00083 


3 


0.75 


-0.94738 ± 0.00053 


6.396 ± 0.059 


0.0103 ±0.0013 


0.09691 ± 0.00098 


3 


0.85 


-0.93846 ± 0.00060 


6.138 ±0.017 


0.0113 ±0.0010 


0.1284 ±0.0015 


10 


0.45 


-4.3575 ± 0.0028 


21.5 ± 1.6 


0.0038 ± 0.0019 


0.0424 ± 0.0020 


10 


0.6 


-4.300 ± 0.015 


19.9 ± 1.4 


0.0044 ± 0.0041 


0.0657 ±0.0017 


10 


0.8 


-4.212 ± 0.034 


20.6 ±3.1 


0.0038 ± 0.0071 


0.1183 ±0.0061 



TABLE II. ISCO configuration state parameters. Results are obtained by minimizing the binding energy per unit rest mass 
Eh for constant rest mass sequences (see Figures ^, ^ and ^ at the highest numerical resolution {rix — 256). Numerical errors 
are calculated by minimizing Et for other numerical configurations (see Table |) and calculating the error as per Eq. Shown 
are tabulated values for the ISCO configurations of the binding energy per unit mass Eb, the total angular momentum of the 
system J, the orbital angular velocity parameter SI, and the maximum rest mass density of the NS {po)^ax various choices 
of BH/NS mass ratio parameter ^tn and total NS rest mass Mq. 





Mo/(Mo)_ 


XA 


Cb 


a 


pV 


1 


0.6 


0.076 ±0.012 


0.87244 ± 0.00049 


2.215 ±0.014 


-0.01288 ±0.00016 


1 


0.75 


0.164 ±0.033 


0.8342 ± 0.0021 


2.113 ±0.059 


-0.01797 ±0.00073 


1 


0.9 


0.2488 ± 0.0098 


0.78917 ± 0.00044 


2.000 ±0.012 


-0.02344 ± 0.00032 


3 


0.75 


0.484 ± 0.037 


0.7886 ± 0.0074 


3.18 ±0.19 


-0.0362 ± 0.0028 


3 


0.85 


0.540 ± 0.021 


0.7605 ± 0.0052 


3.16 ±0.11 


-0.0424 ± 0.0023 


10 


0.45 


0.66 ±0.14 


0.824 ± 0.038 


5.6 ± 1.2 


-0.0274 ± 0.0075 


10 


0.6 


0.73 ±0.12 


0.787 ± 0.068 


5.459 ±0.058 


-0.039 ± 0.022 


10 


0.8 


0.85 ±0.18 


0.76 ±0.13 


9.8 ±3.2 


-0.044 ± 0.050 



TABLE III. ISCO configuration state parameters. Results are obtained by minimizing the binding energy per unit rest 
mass Eb for constant rest mass sequences (see Figures ^ ^ and ^) at the highest numerical resolution {n^ = 256). Numerical 
errors are calculated by minimizing Eb for other numerical configurations (see Table |) and calculating the error as per Eq. [js] . 
Shown are tabulated values for the ISCO configurations of the separation parameter xa, the constant Cs on the right hand side 
of Bernoulli's Eq. ^ the scaling parameter ct, and the momentum of the BH Pg^ (related to the value in the computational 
coordinates {i'} by Pg^ = c^Pbh) f'^'' various choices of BH/NS mass ratio parameter /ib„ and total NS rest mass Mq. 
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FIG. 8. The total angular momentum of the system vs. the 
angular velocity for constant rest mass sequences where the 
ratio of the BH mass to the NS mass, ntn, is 1. Shown are 
sequences whose rest mass are 0.45, 0.6, 0.75, and 0.9 times 
the maximum stable rest mass of a NS for a F = 2 poly- 
tropic equation of state. In units where G = c = fc = 1, 
the maximum rest mass for a F = 2 polytropic star is 
i^o)^ax = 0.179862. Plotted are the results from the highest 
resolution configuration {rix = 256) from Table |, with error 
bars signifying estimates of the numerical error as per Eq. 



The total angular momentum of the system about the 
center of mass (which corresponds to the origin of our 
coordinate system) is defined as ||l^ 

J ADM = ^ / iXj{Kkr,r - K-fkrn)) , (50) 

With our configuration, the z-component reduces to 

xbhP"bh+ J d^x(j)^°pohW^{xvy -yv'') 



J = J 



ADM 



(51) 

Note that the integral in Eq. |5| has support only where 
po ^ 0, namely, inside the NS. 

It is interesting to plot the total angular momentum 
of the system about the center of mass, as was done for 
the binding energy per unit mass in Figures ^ - R In 
Figures || - |l^, we plot the unitless quantity J/-Mj for 
values of mass ratio parameter /if,„ — 1,3, and 10, re- 
spectively. We use the same numerical configurations 
as listed in Table |, and compute the error bars as per 
Eq. In it is argued that, for constant rest mass se- 
quences, the minimum of the angular momentum J and 
the minimum of the binding energy per unit mass Eh 
should coincide, citing a result [ p^ relating parameters 
of stationary solutions to the Einstein equations not con- 
taining black holes, namely, (cIMadm) = ^ (dJADAi)- 
In the study of quasiequilibrium corotating NS/NS bi- 
naries j|] it was shown that the minima of the binding 
energy and angular momentum coincided to "numerical 



accuracy" . It could be asked if the configurations com- 
puted here share this feature. For example, docs the 
minimum of the binding energy per unit mass Eh for 
libn = 3, Mo/{Mo)^^^ = 0.75 (see Figure |) correspond 
to the minimum of the of the angular momentum J for 
the same case (see Figure §)? From Table || we see 
that this minimum corresponds to an orbital angular fre- 
quency of flMo = 0.0103 ± 0.0013. Performing a similar 
calculation for the Mo/{Mo)^^^ = 0.75 case of Figure ||, 
computing numerical errors in the same fashion as for 
Tables || and |m, we obtain a value of orbital angular 
frequency SlA/g = 0.01376 ± 0.00017 at the minimum. 
We see that the orbital angular frequency for the mini- 
mum points do not agree within numerical errors. The 
minimum for the angular momentum occurs at a higher 
orbital angular frequency (thus, at a smaller separation) 
than the minimum of the binding energy per unit mass. 
In fact, for all of our configurations we find it to be true 
that the minimum of the angular momentum for a con- 
stant rest mass sequence occurs at a smaller separation 
than the minimum of the binding energy per unit mass. 

Why is this result different from previous studies 
There are many possible reasons, all of which could be 
contributing to some degree. One difference is that our 
configurations contain a black hole. In this case, the 
result relating stationary solutions to Einstein's equa- 
tions II2I is 



dMA 



1 



DM 



-kh dAn + dJn + ^ dJA 



DM 



(52) 



where kh is the gravitational acceleration on the event 
horizon. Ah is the area of the event horizon, f2jj is the 
angular velocity of the BH, and Jr is the angular mo- 
mentum of the BH; all of these quantities are computed 
on the event horizon of the BH. However, it is certainly 
not clear whether Eq. |5| even applies in our case. It 
was derived assuming a strictly axisymmetric space- 
time (e.g., an axisymmetric BH with matter distributed 
axisymmetrically) , whereas we have a nonspinning black 
hole located away from the center of mass of the configu- 
ration. Another difference is the ansatz used here for the 
extrinsic curvature; Eqs. |l^ and |l^ is slightly different 
than that used in other NS/NS sequence studies 
where the time derivative of the conformal metric is set 
to zero. In order to keep the momentum constraints reg- 
ular at the black hole puncture, we were induced to use 
Eqs. |l| and |l| as the form of the extrinsic curvature, 
noting that this ansatz reduces to the usual form (the 
time derivative of the conformal metric vanishes, Eq. |22| ) 
in the limit as Mbh 0. Finally, we note that the 
present study most likely has a better numerical accu- 
racy, as we have used a resolution that is 2 times finer 
than in It is possible that with both an increase in 
resolution and a better estimate of the numerical error, 
the study in ||] would have detected a difference between 
minima of the binding energy and minima of the total 
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FIG. 9. The total angular momentum of the system vs. the 
angular velocity for constant rest mass sequences where the 
ratio of the BH mass to the NS mass, /itn, is 3. Shown are 
sequences whose rest mass are 0.45, 0.6, 0.75, and 0.85 times 
the maximum stable rest mass of a NS for a F = 2 poly- 
tropic equation of state. In units where G = c = fc = 1, 
the maximum rest mass for a F = 2 polytropic star is 
i^o)^ax = 0.179862. Plotted are the results from the highest 
resolution configuration (rix = 256) from Table |, with error 
bars signifying estimates of the numerical error as per 




FIG. 10. The total angular momentum of the system vs. 
the angular velocity for constant rest mass sequences where 
the ratio of the BH mass to the NS mass, ntn, is 10. Shown 
are sequences whose rest mass are 0.45, 0.6, and 0.8 times 
the maximum stable rest mass of a NS for a F = 2 poly- 
tropic equation of state. In units where G — c = k = 1, 
the maximum rest mass for a F = 2 polytropic star is 
(A'/o)™^^ = 0.179862. Plotted are the results from the highest 
resolution configuration (n^, = 256) from Table |, with error 
bars signifying estimates of the numerical error as per Eq. 0. 



angular momentum. 

VI. CONCLUSION 

In this paper, we have described a method to numer- 
ically calculate general relativistic initial data configu- 
rations corresponding to a binary BH/NS in quasicircu- 
lar orbit. We construct a code to carry out the calcula- 
tion, and carefully validate the code using an independent 
residual evaluator convergence test on each configuration 
solved for this paper. Assuming a F = 2 polytrope, we 
construct sequences of constant rest mass, constant black 
hole mass initial data sets. By minimizing an effective 
binding energy, we find approximate ISCO configurations 
for mass ratios of fibn = Mbh/^ns = 1,3, and 10. A 
technique for monitoring the numerical error is presented 
and used in reporting all numerical results. 

These ISCO configurations are only approximate ones 
in that the constant rest mass, constant BH mass se- 
quences constructed here are not solutions to the full 
Einstein equations (although each set solves the Hamil- 
tonian and momentum constraints). This is true for all 
quasiequilibrium studies ||9|-[l^. Only a full numerical 
general relativistic treatment will be able to quantify the 
errors in quasiequilibrium studies, which may be rela- 
tively large as a binary system approaches the ISCO. 

In this study, we have assumed the NS to be corotating 
about the center of mass of the binary system. However, 
it has been shown that a NS in orbit together with a 



BH would need an unnaturally high viscosity to lock the 
spin and orbital angular velocities. It would therefore 
be more realistic to perform the calculation presented 
here with an irrotational p6|~p8[ NS. This added feature, 
along with a realistic cold equation of state, will be the 
subject of a future study. 
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